FERMILA B-Pub-98/295-A 
tiep-ph/9809489 
September 1998 



Classical inflaton field induced creation of 
superheavy dark matter 

Daniel J. H. Chungfj 

Department of Physics and Enrico Fermi Institute 
The University of Chicago, Chicago, Ilinois 60637, and 
NASA/Fermilab Astrophysics Center 
Fermi National Accelerator Laboratory, Batavia, Illinois 60510 

We calculate analytically and numerically the production of superheavy 
dark matter (X) when it is coupled to the inflaton field within the context 
of a slow-roll m^0 2 /2 inflationary model with coupling g 2 X 2 <j) 2 /2. We find 
that X particles with a mass as large as 1000 Hi, where Hi is the value 
of the Hubble expansion rate at the end of inflation, can be produced in 
sufficient abundance to be cosmologically significant today. This means that 
superheavy dark matter may have a mass of up to 10 -3 Mpj. We also derive 
a simple formula that can be used to estimate particle production as a result 
of a quantum field's interaction with a general class of homogeneous classical 
fields. Finally, we note that the combined effect of the inflaton field and the 
gravitational field on the X field causes the production to be a nonmonotonic 
function of g 2 . 

PACS number(s): 98.80. Cq, 4.62. +v, 95.35. +d 



*Elcctronic mail: djchung@feynmain.physics.lsa.umich.edu 



I. INTRODUCTION 



The rotation curves deduced from observing luminous matter (see for example [|T]) 
indicate dark matter (DM) exists around galaxies. Furthermore, comparison of the 
peculiar velocities of many galaxies with the detailed maps of density contrast suggest 
that Q > 0.3. However, these cannot be all in the form of baryons according to 
the standard scenarios of big bang nucleosynthesis || ||. Structure formation studies 
indicate that relativistic dark matter is unlikely to make up most of the DM |J. These 
evidences suggest the existence of a cosmologically significant abundance of nonbaryonic 
weakly interacting massive particles (WIMPs). Since SUSY models (including string 
inspired ones) generically predict new stable weakly interacting particles, the existence 
of WIMPs is even more likely. 

Despite the fact that the nature of the DM is still unknown, it is usually thought 
that DM particles cannot be too heavy. If the WIMP is a thermal relic, then it was 
once in local thermodynamic equilibrium (LTE) in the early universe, and its present 
abundance is determined by its self-annihilation cross section. As argued by Griest 
and Kamionkowski ||, the self annihilation cross section has an upper bound of ~ 
l/M| from considerations of unitarity, while the temperature at which freeze out occurs 
increases as the cross section is decreased. Hence, the assumption of LTE gives an 
upper bound of about 500 TeV to the mass of the dark matter. The present abundance 
of non-thermal relics (those that never attained LTE) is not determined by their self- 
annihilation cross section because their final abundance is not simply determined by the 
usual freeze out scenario. An example of a non-thermal relic is the axion, and the present 
axion abundance is determined by the dynamics of the phase transition associated with 
symmetry breaking. Non-thermal relics are typically very light, e.g., the axion mass is 
expected to be in the range 10~ 5 to 10~ 2 eV 0. 
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However, nonthermal relics can have masses many orders of magnitude larger than 
the electroweak scale and can evade the unitarity bound of Ref. |J. These nonthermal 
DM particles have been called superheavy dark matter (SDM) in Ref. ||. SDM scenarios 
have been discussed in conjunction with various production mechanisms (see for example 
H and references therein), with the gravitational production mechanism being arguably 



the least fine tuned || [10], ITlj . In this paper, we will explore the idea of jl2|| , which 
is to produce SDM by the same mechanism that is at work in what has been called 
"preheating" scenarios. 

The main ingredient of the preheating scenarios, introduced in the early 1990s, is 
the nonperturbative resonant transfer of energy to particles induced by the coherently 
oscillating inflaton fields. It was realized that this nonperturbative mechanism can be 
much more efficient than the usual perturbative mechanism for certain parameter ranges 
of the theory ]13], [14], [T5J [IB). The basic picture can be seen as follows. Suppose we 
have a scalar field X with a coupling g 2 <f) 2 X 2 where is a homogeneous classical inflaton 
field. The mode equation for X field then can be written in terms of a redefined variable 
Xk = X k a 3 / 2 as 

Xfc(t) + {A + 2qcos{2t)) X k{t) = (1) 

where A depends on the energy of the particle and q depends on the inflaton field 
oscillation amplitude. When A and q are constants, this equation is usually referred to 
as the Mathieu equation which exhibits resonant mode instability for certain values of 
A and q. In an expanding universe, A and q will vary in time, but if they vary slowly 
compared to the frequency of oscillations, the effects of resonance will remain. If the 
mode occupation number for the X particles is large, the number density per mode of 
the X particles will be proportional to |xfc| 2 - If A and q have the appropriate values 
for resonance, Xfe wm grow exponentially in time, and hence the number density will 



attain an exponential enhancement above the usual perturbative decay. This period 
of enhanced rate of energy transfer has been called preheating primarily because the 
particles that are produced during this period have yet to achieve thermal equilibrium. 

This resonant amplification leads to an efficient transfer of energy from inflatons to 
other particles which may have stronger coupling to other particles than the inflaton, 
thereby speeding up the reheating process and leading to a higher reheating temperature 
than in the usual scenario. Another interesting feature is that particles of mass larger 
than the inflaton mass can be produced through this coherent resonant effect. Such a 



process is negligible in a conventional scenario of reheating [17]. This has been exploited 



to construct a baryogenesis scenario [12] in which the baryon number violating bosons 
with masses larger than the inflaton mass are created through the resonance mechanism. 
A natural variation on this idea is to produce SDM by the same resonance mechanism 

HI- 

Interestingly enough, what we find in our work is that in the context of a slow-roll 
inflation with the potential V(<p) = m|0 2 /2 with the inflaton coupling of g 2 <f) 2 X 2 /2, 
the resonance phenomenon is mostly irrelevant to the production of SDM because too 
many particles are produced when the resonance is effective. For the tiny amount of 
energy conversion needed for SDM production (tiny means ~ 10~ 17 of the total energy), 
the coupling g 2 must be small enough (for a fixed Mx) such that the motion of the 
inflaton field only at the transition out of the inflationary phase generates just enough 
nonadiabaticity in the mode frequency to produce SDM. The rest of the oscillations, 
damped by the expansion of the universe, will not contribute significantly to particle 
production as in the resonant case. In other words, the quasi-periodicity necessary for a 
true resonance phenomenon is hardly existent for the case when only an extremely tiny 
fraction of the energy density is converted into SDM. Of course, if the energy scales are 
lowered such that a fair fraction of the energy density can be converted to DM without 



overdosing the universe, this argument may not apply. However, in this paper, we will 
be mostly interested in producing SDM with masses larger than the inflaton mass within 
the context of a large-field inflationary scenario, where this argument will apply. For the 
study of cases in which the resonance starts to become efficient, we refer the reader to 



to Refs. ITq, 15, 13, 12 and references therein. 



The main findings of this work are the following: We find that superheavy dark 
matter with a mass as large as 10 3 Hi, where Hi is the value of the Hubble expansion 
rate at the end of inflation, can be produced in sufficient abundance to be cosmologically 
significant today. Typically, Hi can be as large as 10 13 GeV, which means that the dark 
matter may have masses of the order of the GUT scale. In the process of finding this 
estimate, we derive a simple formula (in the spirit of Ref. |19||), Eq. (fH]), that can be 
used to estimate particle production resulting from a general class of interactions with 
a time varying homogeneous classical field (including the gravitational field). Finally 
we observe that coupling X to the inflaton field can actually decrease the amount of 
SDM produced as a consequence of the inflaton field variation canceling some of the 
nonadiabaticity of the expansion rate responsible for gravitational production of SDM. 

This paper is organized as follows. In the next section we will specify the model and 
the inflationary scenario in which our estimations are carried out. In section III, we derive 
a general formula to estimate particle production when a quantum field interacts weakly 
with a general class of homogeneous classical fields. We then compare our approximations 
to two exact solutions. Section IV follows where we apply the estimation to our model 
described in Section II. We then present numerical results for comparison and a better 
estimation of the maximum cosmologically interesting SDM mass in our model. Finally, 
we conclude with a summary in Section V. 
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II. MODEL 



Two conditions are necessary for the viability of the SDM scenario ||: a) their inter- 
action rate must be sufficiently weak such that local thermodynamic equilibrium (LTE) 
was never obtained and b) the X particles must be cosmologically stable. As we will 
see, because LTE necessitates the reaction rate to be larger than the Hubble expansion 
rate while the reaction rate involves at least an inverse mass squared suppression coming 
from the cross section involved, large mass particles can naturally evade LTE. 

Let us denote px as the energy density of the SDM particles and iix(t e ) as the number 
density of the SDM at time t e when inflation ends. As shown in ||, today's SDM density 
VL X = px(t )/p c (t ) (where p c (t ) = 3H$M P \ 2 /8tt and H = 100 h km sec" 1 Mpc" 1 ) 
can be expressed as 

V To J 3 \M Fl ) M Fl H? y ] 

where Hi is the Hubble expansion rate at the end of inflation, To is the temperature 
today, T-pj^ is the reheating temperature, and fl R h 2 « 4.31 x 10~ 5 is the fraction of 
critical energy density that is in radiation today. Throughout this paper, we will give 
our results in terms of Qxh 2 /S where we have defined 

S = (T RH /(10 9 GeV))(^/(10- 6 M pl )) 2 . (3) 

For a typical reheating temperature of 10 9 GeV, Eq. (||]) implies that the SDM energy 
density today will be Q x h 2 ~ 10 17 (px(^e)/p(^e)) where p(t e ) is the total energy density 
at the end of inflation. It is indeed a very small fraction of the total energy density that 
needs to be extracted to saturate the upper bound on the cosmological mass density. 
Hence, the difficulty of our scenario lies in creating very few particles, if these are to be 
the SDM. 
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Now, consider the nonthermalization condition 

n x (a A \v\)<H (4) 

which allows the evasion of the unitarity upper bound on the mass of DM. Using Eq. 
(□) with n x h 2 < 1 and the fact that for WIMPS, the averaged annihilation cross section 
(ca|i>|) is less than Mj^ 2 (unitarity bound), we can obtain the estimate 

nxMvl) 7xl0- 19 (gg/jfri) 

Hi (T RH /10 9 GeV) (M x /M pl )3 1 ) 

which is the quantity that must be less than one at the end of inflation to avoid ther- 
malization. If Hi « 10~ 6 Mpj and 

^ i jkH_y /s > 10 - a (6) 

\ H, } [lIPGeV ) ' W 

there is no LTE, and the particles density will evolve trivially as was assumed in Eq. ([|). 
Thus, because of the M^ 2 generically coming from the cross section, SDM will generally 
fail to achieve LTE irrespective of the exact value of the weak coupling constant. Note 
that this is a rather conservative estimate since the reheating temperature is likely to 
be larger and the cross section is likely to be smaller. We also remark that because the 
reheating temperature is usually much smaller than the X mass in SDM scenarios, the 
thermal production of the X particles is usually negligible.^] 

For the X particles to serve as DM, they must have a lifetime that is longer than 
the age of the universe and be extremely massive. One possible source of SDM is the 
secluded and the messenger sectors of the gauge mediated SUSY breaking models where 
SUSY can be broken at a large scale (giving rise to large masses) while the secluded and 
the messenger sectors can have accidental symmetries analogous to the baryon number 



1 Since for times larger than t e , the interaction rate continues to be smaller than H , the particles will 
not thcrmalize later either. 
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giving the particles stability []2T, 22]. Other natural possibilities include theories with 
discrete gauge symmetries |3| and string/M theory ||25|| . 

To explore the dynamics which we believe is typical towards the end of inflation, we 
primarily focus on two coupled scalar fields in an expanding universe. The action can 
be written as 



S — S g + Sm 



(7) 



where 



Sm 



d X\ 



d 4 x\/—g 



M pl 2 R 

167T 



— m 



2 j.2 



+ 



g^X^-imjc + ZR + gY)* 



(9) 



We will take £ = 1/6 corresponding to conformal coupling to gravity although our main 
results will be insensitive to this assumption. Neglecting the small effects coming from 
the quantum fluctuations, we take the gravitational field and the inflaton field to be 
purely classical fields. Note that we are neglecting other fields which the inflaton field 
needs to couple to in order to reheat the universe. We have numerically verified that as 
long as the reheating or preheating occurs on a time scale that is greater than 5m7 , our 
main conclusions are insensitive to this assumption since the particle production that 
is mainly of interest to us occurs during this time (as we will later see, nonresonant, 
nonperturbative production is of interest for superheavy dark matter production). 

We will consider a metric of the form ds 2 = dt 2 — a 2 (t)dx. 2 . The resulting equations 
of motion for the homogeneous classical fields are 



d 2 4ir . ■ 2 
t: - — W + m 



2 ±2 



+ 3-0 + m\ 








(10) 
(11) 



where we have neglected the dark matter contribution to the energy density. This is a 
good approximation during the time period of dynamical interest. 

Of course, we do not expect any of our results to be sensitive to the initial conditions, 
since our results depend upon what happens towards the end of inflation and afterwards. 
Our results will mainly depend upon the functional form of the potential for the inflaton 
field. In this paper, we will not study this model dependence but will study what we 
consider to be the typical dynamics of such systems. For the sake of completeness, we 
discuss in the Appendix the initial conditions that we use for our study. 

Now let us consider the X sector. With the canonical conjugate to X as a 3 X and 
canonically quantizing this action, the Heisenberg equation of motion is 

X + 3HX - 4tV 2 X + (M 2 X + g 2 <p 2 )X = (12) 
a 2 

where H = a/a is the Hubble expansion rate. We introduce the Fourier convention 

X = I (27^ (afce * %(t) + a * e ~^W) ( 13 ) 

where we have defined hk = X^a and defined the normalization for the annihilation 
operator as [a fc - , aj-.] = 5^\ki — k 2 ). Imposing the canonical commutation condition, we 
obtain the normalization condition 

h k h* k - h* k h k = -. (14) 
a 

The mode equation satisfied by hk is 



h k + Hhk + 



a \ a 



hk = 0. (15) 



In conformal coordinates defined by ds 2 = a 2 (r])(dr] 2 — <ix 2 ), this mode equation becomes 
hl( V ) + iu 2 k hk = (16) 
where oj k = \Jk 2 + [M\ + g 2 (p 2 (r]))a 2 (r]) and the ' derivative is with respect to conformal 



time. 



Now we need to fix the boundary conditions. Because the particle number can be 
constant only for time translationally invariant systems, the no-particle state (the vac- 
uum state) existing towards the end of inflation can be specified only approximately in 
an expanding universe (see, for example, |27[] and |2S|). One method of systematically 
classifying the various inequivalent approximate vacuum states is through the adiabatic 



vacuum 



27j definition. As will be seen later, we will use an effectively infinite adiabatic 



order vacuum boundary conditions by considering the boundary conditions placed at 
arbitrarily far past and future for the nonsingular spacetime that we will consider. 

If we denote h^ 1 as the mode solution with boundary conditions defined at a future 
time rji and kf? as the mode solution with the boundary conditions defined at a time in 
the past 7/0, we can define the Bogoliubov transformation as h!J} (rj) = akh\^ {rf)+[3kh*^ {rf). 
The number density then is given by 

"M'filwW' (17) 



III. STEEPEST DESCENT METHOD 

This section presents a derivation of the analytic estimation (Eq. (fHP ) that will be 
used in conjunction with the numerical work. A reader interested only in its application 
should skip to the next section. Its direct application to our physical system of main 
interest (presented in Section II) will be given in Section IV. The following analysis is in 
the spirit of Ref. |L9|. 

With the definition 

\/2uJk yZ^k 
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the differential equation 

hl + u 2 k h k = (19) 
is equivalent to 

a' k = -^j- exp (21 J u k dr^j (3 k 

f3' k = ^exp(-2ijw k dri)a k . (20) 

Because ui' k /{2uj k ) vanishes at 77 = ±00 (adiabatic in-out region assumption), ot k and f3 k 
become constants there, assuming no singular behavior occurs there. Expanding a k and 
f3 k in an adiabatic series (in powers of derivatives), and using the boundary condition 
a k(Vp) = 1 an d Pk{.Vp) = (equivalent to an infinite adiabatic order boundary condition 
in the limit r] p — > —00 for our restricted class of spacetime) we obtain 

Pk ~ / rf ^ ex P ("2* fu k (v'W) ■ (21) 

to leading order. Note that this approximation should be good as long as (3 k uj' k / {2ojk) 1 
even when uj' k /{2uj k ) > 1. This is certainly true for the cases to which we wish to apply 
this analysis. Our next objective is to obtain an approximation for this integral. Let us 
write 



u k = sjk* + MlC{r,) (22) 

where all the 77 dependence is contained in C{rj) and the radicand is positive definite for 
all real rj. For example, in our model, 

C(n) = (l + £g2) a^) (23) 

which can be thought of as the square of an "effective" scale factor. We will also assume 
that C(rj) is C°° for real rj and analytically continue uj k to the complex plane. Because 
of the squareroot in the exponent, the poles of the integrand in Eq. (^1|) will also be 
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branch points in the complex rj plane. [] We will choose the branch cuts such that they 
go from the branch points to infinity along a path such that the branch points on the 
lower half plane have the cut going towards —zoo and the branch points on the upper 
half plane have the cut going towards ioo. Furthermore, the cut will be taken along the 
curve where the exponential function has an equal modulus. Transverse to the cut, the 
exponential will fall off rapidly. 

The integral over the real axis can be replaced (using Cauchy's theorem) with the 
integral over an appropriately deformed contour in the lower or upper half plane (we will 
soon see that our phase convention is such that we are really concerned with the poles on 
the lower half plane as shown in Fig. |ip. The main contribution from the integral over 
the deformed contour will come from near the branch points and possibly the end points. 
The branch points will be distributed symmetrically with respect to reflection across the 
real axes because of the Schwartz reflection principle. The end point contribution will 
be of the order of u' k /uJk. However, in our restricted class of spacetimes which admits an 
infinite adiabatic order vacua, we can make this contribution arbitrarily small by taking 
the endpoints further out. (We comment further on this effect later.) Hence in general, 
Eq. (plD can be approximated as a coherent sum of steepest descent integrals around 
each of the branch points. 

Let us look at the contribution from the jth branch point denoted as fjj. Near this 
branch point, the integral in the exponent of Eq. (^T|) can be expanded as 

fujMdv = r + ^V^W 372 + - ( 24 ) 

Jr/p Jr/p 6 

where 8 = rj — fjj and we kept the leading term in the 5 expansion of C (we will assume 
that C does not vanish here). Expanding in a similar way oj' k /uJk in Eq. (^1|), we can 
2 We assume that the integral of the squareroot in the exponent will introduce no other branch points. 
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write the contribution to (3k from this branch point as 

Uj=Vjexp(~2i / u)k(rj)drj) (25) 

Jr] p 

where we defined Vj as 

vj = \ J Cj f eM^M x ^W)5^) (26) 
and Uj was introduced to rewrite (3k as 

&«E^- (27) 

j 

Here the contour Cj near ^ is the steepest descent contour. 

Let us determine the steepest descent path near fj. If we denote 9 2 = arg(C"(?7)) and 
denote 6 to be the argument of 5 along the steepest descent contour, the restriction on 

e is 

, = *(*-i>-ft (28) 

o 

where n is an integer. Now, if we let the branch cut go along arg(5) = a such that 
arg(<5) G [a,a + 2n) on the lower half plane and arg(<5) G (a — 2n, a] on the upper half 



plane, upon choosing y C'(?7) > when 77 is real (consistent with positive frequency mode 
definition), we can place the restriction 9 2 G [— a,2n — a) in the upper half plane and 
62 G (— 27T — a, —a] in the lower half plane. As mentioned above, we choose the branch 
cuts to go towards ±200 by restricting a G (0, ir) on the upper half plane and a G (— 7r, 0) 
on the lower half plane. Finally, restricting the branch cut to be on an equal modulus 
curve, we obtain the relationship 

-e 2 ± 2vr 

a = (29) 

where the positive sign corresponds to the upper half plane and the negative sign corre- 
sponds to the lower half plane. 
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Combining Eqs. ( p8|) and (^), we find two possible integers n — 0, 1 for the lower half 
plane, while we find only one possible integer n = for the upper half plane. Since the 
steepest descent approximation requires an incoming direction and an outgoing direction 
for each of our branch points, we can only use the lower half plane branch points. To 
summarize, the necessity of making branch cuts with the appropriate shape and the 
information regarding the derivative of C at the branch point determines the usable 
steepest descent paths, which are restricted to the lower half plane (lower and not upper 
because of the sign convention on the squareroot). This, as we will see below, corresponds 
to an exponentially damped result. 

Let us return to the evaluation of Vj in Eq. (|26|) . From Eq. (p8|) and Eq. (p9|) , we see 
that the branch cut bisects the angle made by the steepest descent contour which makes a 
fixed (acute) angle of 2tt/3 (or 4n/3 in coordinate angle), and hence the steepest descent 
path is well defined independently of 6 2 . Hence, we can easily evaluate the integral after 



making a change of variables u = 5 3 ^ 2 (—AimJC'(fjj)/3). The result is 



Vj = f (30) 



Uk{v)dv (31) 



Now we need to evaluate 

m 

J ri P 

to complete our evaluation of Uj. Since we do not in general know the function C{rj) 
on the complex plane but know it and its derivatives on the real axis (numerically for 
the model presented in section II), let us find an estimation utilizing that uo k is analytic 
around the real axis. Let fjj = r + /i where r is purely real and /i is purely imaginary. 
Now, we split the integral into 
nj 

/ u k (r])dr) = + Jj (32) 
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where 

$j = [ r Mv)dr] (33) 

J tip 

Jj = J u k (r))dr) (34) 
where $j is purely real. To evaluate Jj, we expand in Taylor series, 

Jj « Uk(r)» + ^(r)M 2 /2 + <(r)/i 3 /6 + - (35) 

where one should note that all the even terms are real. Because we will mainly be 
interested in one pole domination case, we will only calculate the imaginary value of Jj. 
Thus, as long as 

KM I « 16/V | (36) 

we can truncate the Jj after the first term. We shall check the self consistency later. 

To approximate fjj, we Taylor expand the left hand side of C — —{k/M x ) 2 about r 
to obtain 

f*C m (r) + C(r) = (37) 
o 

,2 , ,2/ 



T c " w = ~mT (38) 



where C"(r) > since /i is purely imaginary. The truncation of the Taylor expansion 
should be justified as long as 



C"'(r) 



" CT « 1 (39) 



lQM 2 x C"{r 



\CW(r)\ 



ui 



2 



< 1 (40) 



C"(r) 2 6M X 

for Eqs. (J37f) and (^) respectively where the superscript indicates the order of the 
derivative.^ 

3 Note that even if Eq. (|9|) is not satisfied, all the equations derived may still be valid as long as r 
can be estimated another way. 
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Assuming that the contribution from one of the poles dominates in Eq. (|2T|), we can 
approximate \/3k\ 2 ~ |tA| 2 - Using Eq. ( p8|) and choosing the solution on the lower half 
plane (as we have justified above), we obtain 

(k/a eS (r)) 2 M x 



\(3 k \ 2 w exp 



M x ^fHl s (r) + R eS (r)/6 jH^r) + i? eff (r)/6 
where we have dropped the factor of (7r/3) 2 , we have defined the effective scale factor 



(41) 



a eff = VC, (42) 

and -ff e ff (r) and i? e g(r) correspond to Hubble expansion rate and Ricci scalar (respec- 
tively) for the metric ds 2 = CL^(r])(dr] 2 — c?x 2 ). All time dependent quantities in Eq. ( f4T|) 
are evaluated at rj = r given by Eq. ([IT]). Explicitly, the radicand of the exponent in 
Eq. (f4j is simply H* s (r) + R eS (r)/6 = C"(r)/(2C 2 (r)) since 6a^ ff /a3 ff = R qS (r) and 
a eff/ a eff = ^eff^ r )' R- ewr iti n g the condition Eq. ( PS| ) in this new notation as (in the 
k — > limit) 

6H 2 eS (r) + R eS (r) » R eS (r)/6 (43) 

we see that this condition is almost always satisfied as long as the Eq. (|40|) is satisfied.0 
The exponentially suppressed behavior is not exact for arbitrarily large masses since 
the endpoint contribution will eventually become larger. An asymptotic analysis of the 
problem in the appendix of || gives the general power law mass dependence of the 
endpoint contribution. As discussed there, this contribution can become particularly 
significant when C is not analytic. The present analysis, however, shows that owing 
to the fact the endpoints lie far away from the branch points (which have stationary 
phases), in the case that C is analytic, the contribution from the endpoints can be 
made arbitrarily small as long as the spacetime admits an infinite adiabatic order in-out 
regions. 

4 Eq. (^) must be satisfied since we are using Eq. (|3|) to get the value of fi, 
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Let us restate the main point of this somewhat technical section. Given a homoge- 
neous classical field coupling that gives rise to mode equations of the form Eq. ([H]) with 
Uk given by Eq. (p2[) and given an in-out infinite adiabatic order vacua, the Bogoliubov 
coefficient giving the particle density per mode is given by Eq. ([II]) and the number 
density is given byf] 



u\ Q eff( r ) ( - AM x 

n x(t) ~ o/o ^ta exp 1 



M 



x 



3/2 



(44) 



where a e g(r), i? e g(r), and if e g(r) are effective scale factor, Ricci scalar, and Hubble 
expansion rate defined by Eq. (E3) and the statement following it. The "effective" quan- 



tities are all evaluated at a value of r satisfying Eqs. (]3~7f ) and (|38|). As will be discussed 
towards the end of next section, r is roughly the point at which C varies nonadiabatically, 
or a bit more precisely, where C"(r)/C 2 (r) is a maximum. For couplings of the form 
Q(4>)X 2 /2 + £RX 2 /2 where Q is analytic and R is the Ricci scalar of the spacetime, we 
can write 

C{r,) = a 2 ( V )(l + Q{4>)/M 2 X + (£ - 1/6)J2/M£) (45) 
in a spacetime with the metric ds 2 = a 2 (r))(dr] 2 — c?x 2 ). 

IV. COMPARISON WITH EXACT RESULTS 

Let us first apply Eq. (fHP to a couple of exactly soluble cases. Consider the case 



given by Ref. [30] when there is no inflaton coupling to the X field and the spacetime is 



5 Although r in general depends on k (for an explicit example, see Eq. (|5q ) ), because most of the 
contribution to the number density comes from (k/Mx) 2 <C 1 where the k dependence can be neglected, 
we shall neglect k when doing the momentum integral. 
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given by ds 2 = C(r])(di] 2 — dx 2 ) where 



2^2 



C{rj) = cf + cfy] 



where C\ and c 2 are real constants. The exact number density per mode is 



|/3 fc | 2 = exp 



-7T 



+ 



M x c 2 



(46) 



(47) 



M x c 2 c 2 

which has been noted |K| as the spectrum of a nonrelativistic gas of particles with mo- 
mentum k/y/C having chemical potential of —M x c\/(2C) and a temperature of c 2 /(27tC) 
In deriving Eq. (|4l|), we have made two separate approximations. One is the usage 
of the steepest descent method and the other is the estimation of the integral Eq. (|3i"D . 
To test the goodness of each of these approximations, since we know C exactly here, we 
will first consider the steepest descent method with the exact branch point. Let's start 
from Eqs. fl25|) and One branch point is on the lower half plane at 



7] = 7] 



\ 



c\ + k 2 /M 2 x 



r 2 



(48) 



and another, its complex conjugate, is on the upper half plane. The quantity Vj has been 
already evaluated in general, corresponding to a contour integral around the branch point 
on the lower half plane, and is given by Eq. (p0|). The integral Eq. (|3l| ) can be easily 
calculated with fjj — fj to give 



/ u>k(r])dr] = (—in/ 2 + real phase). 

Jrjp 



Putting this in Eq. (|2q) , we have 



71 



exp 



-it 



+ 



M x c 2 



_M x c 2 c 2 

whose exponent we see matches the exact result of Eq. 

Now, let's apply Eq. fl41|). Eq. (p7|) gives r = 0. Hence, we find 



(49) 



(50) 



exp 



k 2 M x c\ 



M x c 2 



C2 



(51) 
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which gives only an approximately correct exponent. However, note that the functional 
form of the mass dependence of the exponent is exact. Hence, although the steepest 
descent approximation seems to give accurate results for the exponent, as exemplified by 
Eq. (|5DD, the Taylor expansion method used to estimate the integral Eq. ( |3"TD leads to a 
numerical value of the exponent that is only roughly correct as seen in Eq. (fTTD . Still we 
see the functional behavior of the mass dependence is accurate. 

Let's consider another exactly soluble case |3l| (see also f28|) specified by 



C(r}) = A + Btanh(fyq) 



(52) 



where A, B, and p are positive constants with A > B, This results in particle density 
per mode of 



sinh" 


/ w 


^Jk 2 + M 2 X {A + B) - 


\Jk 2 + M X (A - B) 


} 


sinh 


^k 2 + M 2 x (A-B)_ 


sinh 


^k 2 + M 2 x (A + B) 



(53) 



When the exponential cut off starts to become effective for large masses, the behavior is 
W * exp (-^LVA^b) (54) 



P 



where we have effectively set k = 0. Let us compare this with the steepest descent 
approximation. In this example, it is easier to find the branch point exactly. However, 
since we are interested in using the Taylor expansion estimation, let us solve Eqs. (|37D 
and (|38|) which are valid if the conditions Eqs. (|39|) and (flOl) can be satisfied. These 
conditions are respectively 

-20Bp 2 (l-t 2 r )t r 



4p 2 (2 - 17t 2 + 30tj - 15t 6 r ) 



;i-^)(i-3t 2 ) 

Ap 2 t r (2 - 5t 2 + 3tf) 



(l-t 2 )t r 



(k/M x ) 2 + A + Bt r 
-12Bp 2 (l-t 2 )t r 
(k/M x ) 2 + A + Bt r 



(55) 
(56) 



where we have defined t r = tanh(pr). Here, we know that t r < because 



C"{r) = -2Bp 2 {\ 



(57) 
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must be positive as we explained before. Hence, we only need to consider t r G (—1,0] 
and find that our conditions imply A ~ B. Solving Eqs. ( p?|) and ( |38| ) for t r , we find 



—B - ^/ffl + 3(A + {k/Mxff 

3(A+(k/M x y) (58j 

which then implies 

^ 4 + 0(A- B , (59) 



2C 2 (r) 9(A - E) 25 
where we have expanded in powers of A — B and set k = 0. Therefore, in the large mass 
limit, using again Eq. (f41~l), we have 

|/3 fc | 2 « exp L^VA^b) (60) 



which is in reasonable agreement with Eq. (|54|). The lesson that we have learned here is 
that as long as one accounts for the regime in which the approximations are valid, the 
parameter dependence of the damping exponent is accurately calculated by our method. 

Finally, we note that an intuitive meaning can be attached to the value of r. The 
quantity H^(r]) + i? e g(^)/6 has a maximum at an ^ = rj* satisfying the equation 

C'C" + CC" = 3C'C" (61) 

whichP] is very similar to the equation determining r 

{k/M x ) 2 C" + CC" = 3C'C". (62) 

Since {k/ Mx) 2 C" ~ C'C" (as can be seen by looking at Eqs. ( |BTD and ([ill) ) for important 
values of k and since r/* is a stationary point, we expect r to be very close to r/*. Therefore, 
r in general will roughly correspond to the most "nonadiabatic" point, i.e. the point at 
which H^(t]) + R g q(t])/6 is a maximum. 

6 Taking into account the validity condition for the Taylor expansion, it is easy to see that 77* indeed 
corresponds to a maximum. 
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For example, let us consider the model specified by Eq. (0). If one takes the branch 
point determining equation to be Eq. (|6T|), then the resulting value of the number density 
per mode is 

W « ex P (-^^VA^j (63) 
which is in reasonable agreement with Eqs. (|60"D and fl54j). 



V. ANALYTIC ESTIMATION 



Let us first consider the work of Ref. [T8|] to see what part of the parameter space we 



are interested in. Ref. |L8| estimated the particle production in our type of system by 
neglecting any contribution from the decaying mode and using the Mathieu instability 
band plot to estimate the "average" exponent of the mode growth. (Recall that in an 
expanding universe, the mode equation is not exactly the Mathieu equation, since the 
parameters of what would be the Mathieu equation is time dependent in an expanding 
universe.) This method gives a good idea of when the particle production will be efficient 
through the resonance phenomenon. 

Let's start with the number density per mode written as 

\Pk\ 2 = ^ (W + hx\ - i-Xk\ 2 ) - 1/2 (64) 



where Xk is the solution to the conformally coupled mode equation (jl5f) with Xk = \^h k 
satisfying the Oth order adiabatic vacuum boundary condition in the past and Q k — 



(k/a) 2 + M\ + g 2 (f) 2 . Approximate the relevant solution as 
exp(/ fim^dt - i J w k dt) 



where w k = y^l + (l/4)i/ 2 — (1/2)H and \i is the Mathieu characteristic exponent 
when the mode equation is written in the form of Eq. (j^) (see below). The form of the 
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approximation is roughly equal to the lowest order adiabatic approximation and taking 
only the growing mode is justified by the fact that the decaying solution will be small in 
comparison. Now, approximating Wk ~ and assuming (/xm^) 2 / (2S7 2 ,) <C 1, we obtain 
the approximation 

|/3 fc | 2 ~i[exp(2| / um^)-l] (66) 
which is the starting point in ||18|1 f\ 



The mode equation flT5|) with Xk = when written in the form of the Mathieu 

equation Eq. ([I]) gives 



A » ( _L V + ( Ml V + 2q + 2/(3m,t) 2 



(67) 



where we have chosen the initial time to be at t = i = l/ m </> and $(to) represents 
the amplitude of the inflaton field oscillations at the initial time when the inflaton field 



oscillation amplitude starts to decay like 1/t in a pressureless universe. Ref. |18| then 

considers the trajectory of these parameters as a function of time and estimates the 

exponent integral in Eq. (|66|). They find that the efficient preheating ceases at around 

Mx = lOOm^, for g < 1. After that, the incoherent decay process dominates the energy 

release of the inflaton. 

In our case, we are more interested in the regime in which the resonance phenomenon 

is extremely ineffective, or virtually non-existent, because only an extremely tiny density 

of very massive particles is allowed for the dark matter scenario. Define the parameter 

b = Mx/Hi where Hi is the value of the Hubble parameter at the end of inflation. 

Consider the case when 6 = 7 and g = 0, for which Vtxh 2 / S = 0.128. When the 

7 They had used minimal coupling to gravity, but as we can see that is of little importance as far as 
the approximation Eq. ([36]) is concerned. 
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interaction is turned on to say gMpi/ifj = 45, we find that there is a factor of 100 
increase in the particle density produced to Qxh 2 /S = 12. As illustrated in Fig. |2|, 
we see that for this value of gMp^/Hi, b, and values of k that dominate the integral 
for the mode sum, the Mathieu parameter trajectories never cross any of the instability 
bands. Thus, the analysis of Ref. [18] is only partially applicable to our case, mostly 



because we are only requiring a very small fraction of the inflaton energy density to 

turn into X particles. This means that the X particle masses that can be produced in 

interesting quantities should be significantly greater than lOOm^, for coupling constants 

still less than 1. Furthermore, since significant dark matter production occurred without 

the trajectories ever crossing any of the instability bands, we can expect sufficient dark 

matter to be generated without any resonance effects. 

Now let us apply the estimation of Eq. ( [41"| ) to the system presented in Section II. 

Let us consider the largest possible perturbative coupling of around gM-p\/Hi = 10 6 

which is what would give the largest possible mass for the dark matter produced.^ We 

first look for the solutions to Eqs. ([J7]) and fl3"B|). There are many solutions to these 

equations, as we expect. However, only one of these solutions will be relevant for the one 

pole domination approximation. To see which solution will be relevant, we plot in Fig. 

H] the absolute value of the exponent in Eq. ([14]) which we will denote ctS 7~, clS db function 

of time values that are near the actual solution to Eqs. (|37D and (|38|). For definiteness, 

we will fix b at 1080 although the one pole domination approximation will be reasonable 

for all other masses within the range of our interest. It is clear from these plots that 

the pole having a real part corresponding to t m 109 will dominate in the sum Eq. (f27|), 

thereby justifying the one pole approximation reflected in Eq. ([11]). In Fig. |3|, we have 

also plotted the case of gMpi/ifj = 10 3 and b = 30 to indicate the generic nature of this 

8 We are assuming here that something like SUSY is protecting the inflaton potential from large 
radiative corrections that can spoil inflation. 
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one pole domination. 

The solutions to Eq. ( ftTD and Eq. (|38| ) that is of interest to us are not sensitive to 
the value of b = Mx/Hi as long gMpi/ 'Hi/b ^> 1 which turns out to be within our range 
of interest mainly because < M-p\ and \<p/Hi\ < M-p\ at the end of inflation as can 
be seen by looking at the exponent of Eq. (|4l|) . For example, with k = 0, the value of 
t(r) that will be of interest to us is 108.90 (in units of 1/(12.2^)) for both b = 900 and 
b = 1800. Only for values of b as large as 10 4 (which we will not be of much interest to 
us), does the value of t(r) deviate to 108.89. The dependence on k, which can be seen 
directly in Eqs. (|37|) and fl38|) will also be negligible because the most of the SDM that 
are produced will be nonrelativistic. 

Hence, with t{r) = 108.9 we use Eq. (|TD to find that Vt x h 2 / S = 1 at b « 1450 and 
Qxh 2 /S = 0.01 at b ~ 1550. Hence, taking ifj « itl^/2, the maximum possible mass of 
SDM for which its abundance will be cosmologically significant is an order of magnitude 
above the maximum mass that can accommodate efficient preheating. In the next section 
we will give a bit better estimate by taking into account the numerical calculation of the 
particle production. 

Now, let us check the validity of truncating the Taylor expansion in approximating 
the location of the poles. The quantities on the left hand side of Eq. ( flOf ) and Eq. 
([39]) are plotted in Fig. From the figure it is clear that for the gM-p^/Hi = 10 6 the 
approximation is well justified. However, the right hand panel of the figure shows that 
the approximation is only marginally adequate for gM-p^/Hi = 10 3 . 

VI. NUMERICAL ESTIMATION 

Because \f3k\ 2 required to have Qxh 2 / S ~ 1 is extremely small and one must in general 
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compare oscillatory functions to obtain it, it is difficult to calculate it numerically. For 

Qxh 2 /S ~ f, one requires an accurate calculation of 
f CT 5 

^ ~ m 

where k/(a,iHi), the average momentum component, is typically around 0.36. Hence, 
for b ~ 10 6 , even with appropriate scaling, the calculation is numerically delicate. Fur- 
thermore, since many momentum components must be calculated for the integration of 
the spectrum to obtain the number density, an accurate calculation is time consuming, 
at least within a straight forward framework of calculation. The results presented in 
this section comes from a Runge-Kutta solution to a system of equations including Eqs. 
(|T5|), (]T0|) , and (|TTD, all appropriately scaled both in the independent and the dependent 
variables. 

A sample of the numerical results are presented in Figs. |5| and []. Motivated by our 
experience with the comparison with exactly soluble cases in section IV, where we have 
seen that the Taylor expansion approximation used gets the exponent correct only up 
to a constant multiplicative factor close to 1, we introduce a correction factor / in the 
exponent of Eq. (f|l]) as 



|&| 2 «exp -4/ 



(k/a eS (r)) 2 + M x 



(69) 



Mx V^dffW + R eS( r )/ Q V^ffW + ^eff( r )/ 6 . 
and integrate this with respect to k to get Qxh 2 /S that is plotted with the solid curve in 

Figs. [5]and|6|. It is clear from the fact that / is about 1 in both the gM-p\/Hi = 10 6 and 

gM-py/ 'Hi = 10 3 cases that our analytic approximation is a reasonable one. Furthermore, 

as exemplified in Fig. [F], the fit to any particular spectrum by adjusting only / gives 

nearly an exact fit with / between 0.9 and 1.0. From Fig. [5], Mx ~ 1700fZj ~ 10 3 m ( f > 

is the maximum possible value of superheavy dark matter that can be produced in 

cosmologically significant abundance in this inflationary scenario unless Hi > 10" 6 M pl 

or T rh > 10 9 GeV. 
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Before concluding, we would like to note a phenomenon that occurs when the nona- 
diabaticity of the mode frequency change caused by the inflaton coupling is comparable 
to the nonadiabaticity caused by the gravitational coupling. If we crudely approximate 
the exponent of Eq. (^Tj) as 



-4 



(k/a eS (r)Y M x 
M x H eS (r) H eS (r) 



(70) 

which is qualitatively reasonable in our scenario, then we can find using crude estimates 
made in the Appendix (i.e. Hi pa <\> ~ Mpj/Vl27r, ari d <p pa — Mpp^/(2\/37r)) 

that Vtxh 2 / S vanishes at gM-p\/Hi pa 46. Of course, in reality the nonadiabaticity is 
not perfectly canceled, and we expect only a dip in the particle density as gM-p\/Hi 
is increased from 0. Indeed, this cancellation of nonadiabaticity has been observed 
numerically as illustrated by Fig. |8|. Hence, in general, for small positive g 2 couplings 
to the inflaton field, the particle production is not a monotonic function of the coupling 
constant because of the presence of the classical gravitational field. 



VII. SUMMARY 



In this article, we have considered the production of superheavy particles X that 
are coupled to a homogeneous classical inflaton field <fi. We have found that within the 
context of a reasonable V{4>) = m^(j) 2 /2 slow-roll inflationary scenario with the coupling 
g 2 X 2 <j) 2 /2, the parametric resonance phenomenon tends to overproduce the number of 
dark matter particles. Only when the resonance phenomenon completely shuts off is the 
number density small enough to be consistent with the constraint Q x h 2 < 1. For the 
superheavy particles to be cosmologically significant today (Qx ~ 1), its mass can be as 
large as about lOOOifj ~ lOOOm^, as indicated by Fig. |5]. Taking into account the COBE 
normalization of the curvature perturbation power spectrum, i.e. 2Pj/ 2 /5 ~ 1.91 x 10~ 5 
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at a scale of k ~ 7.5H (see for example and using the approximation 



2 1/a = 8^0F V3/ 2 

5 R V75M p f V 1 J 

where the inflaton potential is evaluated at the time of horizon exit (about 50 e-folds 
before the end of inflation), we can estimate ~ 10 13 GeV. This means that the mass 
of the dark matter particles can be of the order of the GUT scale. 

In the process of making the superheavy dark matter mass range estimate, we have 
derived a simple general formula Eq. (|4lD giving an estimate of the particle produc- 
tion due to interactions with classical fields in the limit that there is only one, dom- 
inant nonadiabatic time period during which particle production "occurs." The time 
dependent quantities in Eq. fl4T|) can be approximately evaluated at the point at which 
C"(r)/C 2 (r) (primes refer to conformal time derivatives) is a maximum, where for ex- 
ample, for couplings of the form Q((p)X 2 /2 + £RX 2 /2, C is given by Eq. (|45|). This 
result is applicable to almost any case of time varying homogeneous classical scalar field 
interacting with a quantum field as long as the number of particles produced is small. 

Finally, we pointed out a phenomenon in which the number of X particles produced 
actually decreases as the coupling to the inflaton field <fi increases (Fig. ||). This is a 
simple consequence of the fact that the nonadiabatic variation of the inflaton field is 
canceled out by the gravitational effect of nonadiabatic change in the scale factor. 

As far as the observability of SDM is concerned, the prospects seem no better than 
in the case of electroweak scale WIMPs, unless the WIMPs are strongly interacting or 
charged. Previously, charged p3| , |34[| or strongly interacting dark matter []35| has been 



ruled out with a combination of arguments coming from the unitarity bound |6j and 
experimental observations. Since the SDM evades the usual unitarity bound, charged or 
strongly interacting dark matter may still be viable in this scenario. 

If the SDM decays via an electromagnetic or a hadronic decay, its decay products 
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may change the spectrum of the diffuse gamma ray background [[JIJ . Hence, the diffuse 
background photon measurements of EGRET and COMPTEL give strong constraints to 



such decaying particles for a fairly large range of life times. Note that Refs. [[U], [37], |38 
also suggest that if the SDM decays, its decay products may manifest in the form of 
ultrahigh energy cosmic rays. Note that in these decay scenarios, the life time must be 



in general much longer than the age of the universe. For example, according to Ref. |38 
if the cosmic ray energy spectrum above 10 11 GeV is produced by the decay of SDM, the 
lifetime of the SDM must be around 10 12 £x in units of the age of the universe where £x 
is the fraction of the cold dark matter in SDM. 

In general, we do not expect the direct WIMP detection experiments to be sen- 
sitive to SDM because of their low abundance. The detection rate which goes like 
R ~ po<jv I \M x m N ) where p is the matter density of the halo, a is the elastic scattering 
cross section, v is the virial speed, and is the mass of the nucleon. Hence, unless if 
the particles are strongly interacting in which case a ~ lOmb |39]], the WIMP detectors 
will not have sufficient event rates to measure these particles. 

The indirect method of dark matter detection (detecting the energetic neutrinos 
produced by annihilation of dark matter captured through elastic collisions ||40|| ) will 
also have difficulty in the SDM scenario. The neutrino detection rate in general depends 
upon the SDM's capture rate in the sun through elastic collisions, its annihilation rate, 
and the neutrino cross section for the production of leptons in the rock or the detector. 
Since SDM mass will be much greater than that of the elastic scatterer, it will lose very 
little of its momentum per elastic collision (fractionally m^/Mx where mjy is the mass 
of the nucleon). Hence, in addition to the small number density (0.4 GeV/cm 3 /Mx) 
to begin with, the capture probability through elastic collisions will be negligible. The 
annihilation rate will also be suppressed (unitarity bound ~ 1/M X ) even if one assumes 
maximal branching fraction to the neutrino producing channels. However, the the cross 
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section for the production of leptons in the rock or the detector will be significantly 



enhanced. Still, because the cross section will only grow like \[Mx for Mx much greater 
than the mass of W (assuming that the neutrinos will have energies that scale like Mx), 
the enhancement will not be sufficient to overcome the suppressions. Indeed, even if one 
neglects the neutrino absorption rate in the sun, if there is no significant accretion of 
SDM in the sun, one can easily show that the detection rate will be much smaller than 
the current detector sensitivity |4(| of 10~ 2 m _2 yr _1 . 

Because the dark matter searches have focused mainly on those particles with masses 
that are less than about 100 TeV, the observational consequences for SDM have been 
relatively unexplored. Since as shown in Refs. || |J and in this paper, production mech- 



anisms exist for such particles, and since Refs. |2T], gj, |25| has shown such particles exist 
in extensions to the standard model, a more careful study of observational consequences 
of SDM scenarios may be worthwhile. 
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APPENDIX 



In this appendix we remind the reader of the boundary conditions for the classical 
fields (i.e. Eq. (|TU]) and Eq. (JTTj) ) which is solved numerically. As we mentioned earlier, we 
will choose the initial conditions as to set up a slow-roll, large-field inflationary epoch, at 
the end of which the particle creation will occur. To gain intuition for the needed initial 
conditions, consider the slow-roll scenario as presented in |p6 . From the Friedmann 
equation, 



H 



* " m (Al) 



\ 3M pl M pl 

where 2 <C V{4>) has been assumed (slow-roll scenario) and V(<p) is the inflaton po- 
tential. Since in general <ft ^ 3H<p is required for sufficient inflation to occur, we can 
write 

3if0 « (A2) 



Combining this with Eq. (|A1| ) gives 

/ _8vr /•*(*) V(d>) \ 
a « a(t p ) exp — - / (A3) 



M p f 7^(t p ) 1/^(0) 

where t p (p stands for past) is the time at the beginning of inflation. Hence, for potentials 
of the form constant x <ft n , the number of e-folds of inflation is determined solely by the 
initial and the final values of 0. 

The end of the slow-roll scenario is obtained by determining when the potential energy 
of the inflaton becomes comparable to the kinetic energy. Eq. (|A1[) and Eq. (|A2|) combine 
to give (in addition to Eq. (|A3|) ) 
• _ M m V 4 



(A4) 
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which when combined with the approximate condition for the end of inflation (namely 
4> 2 /2 « V(4>)), we have 

vm e )) „ m p1 (A5) 



where t e is the time at the end of inflation. 

Therefore, for potentials of the form constant x <fr n , fixing the number of e-folds fixes 
the initial value. For V{4>) = rrAcf) 2 /2 potential, we thus have 



0(t e )^M pl / v / 127 (A6) 

and 

M pl 



#*p) w-^V^ + Ve (A7) 



where Af e is the number of inflationary e-folds. For our typical runs, we use N e « 64. 
Then Eq. (|A4j) (equivalent to assuming = 0) gives the initial condition for <fi. 
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Figure 1: A schematic sketch of the analytic structure of Eq. on the complex rj 
plane is shown. The crosses represent branch points and the lines emanating from them 
branch cuts. Shown also is a schematic sketch of the appropriately deformed contour for 
the steepest descent approximation on the lower half plane. 
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Figure 2: Shown is the Mathieu instability chart. For the values of A and q that fall in 
the shaded region, the solution to the Mathieu equation is unstable. The line connecting 
the origin to the upper right corner of the plot is A = 2q. All A and q trajectories must 
lie above this line by their definition in Eq. (|6?D. The three short line segments near 
(q,A) = (0, 11) represents three representative trajectories each having different values 
of k (for those k values that dominate the number density integral) for gM-p\/Hi = 45 
and 6 = 7. Although none of the trajectories cross any of the instability bands, there is 
an enhancement by a factor 100 in the particle density production. 
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Figure 3: Here we plot the number density suppressing exponent T=Ab/(H^(t) + 
i? e g(t)/6) 1//2 as a function of t (in units of ~ 1/(12. 2ifj)) near the values corresponding 
to the conformal time r which solves Eqs. ( |3"7| ) and ([^). The solid curve corresponds 
to the case b = 1080 and gM-p-^/Hi = 10 6 . The dotted curve corresponds to b = 30 and 
g = 10 3 . The actual values of r solving Eqs. (|37]) and (|38"D with these parameters lie very 
close to the minimum of these curves (there exists only one solution r per curve shown). 
Exact numerical value examples are given in the text. The t(r) near 109 and near 126 
gives the smallest two r values for each parameter set. Clearly the pole having a real 
part r with t{r) near 109 dominates. 
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Figure 4: We plot as a function of b = Mx/Hi the quantities that must be less than 
unity for the Taylor series truncation to be valid. The solid curves correspond to the 
left hand side of Eq. ([5]), and the dotted curves correspond to the left hand side of 
Eq. (|39|) . The lower of the given curve type corresponds to k/(a,iHi) = 0. The upper 
of the given curve type corresponds to fc/(ajifj) = 1000 for the gM-p\/Hi = 10 6 case 
and k/(a,iHi) = 30 for the gM-p\/Hi = 10 3 case. These k values indicate the momentum 
range over which the particles are produced. The expansion truncation is clearly justified 
in the gM-p\/Hi = 10 6 case while it is marginally adequate in the gM-p\/Hi = 10 3 case. 
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Figure 5: The solid curve in this plot of Vt x h 2 / S versus b for gM-p\/ Hi = 10 6 shows the 
numerical results. The correction factor (explained in the text) used for the dotted curve 
is 0.932 and for the dashed curve is 0.903. Hence, unless if Hi > 10 _6 Mpj or Tj^g > 10 9 
GeV, cosmologically significant dark matter production in this scenario does not occur 
for masses above about 1700/7^, which can be in the GUT scale. 
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Figure 6: The solid curve in this plot of Qxh 2 /S versus b for gM-p\/Hi = 10 3 shows the 
numerical results. The correction factor (explained in the text) used for the dotted curve 
is 0.98 and for the dashed curve is 0.963, both of which are quite close to the values used 
in Fig. |5|. Hence, we see that our analytic approximation seems to be in agreement with 
the numerical results. 
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Figure 7: For the parameters gM-p\/Hi = 10 6 and b = 1080, the particle density per 
mode is plotted as a function of the wave vector scaled by the scale factor dj and the 
Hubble expansion rate Hi at the end of inflation. The solid curve corresponds to the the 
brute force numerical results. The dashed curve corresponds to the spectrum obtained 
using the analytic approximation Eq. (|69|) with / = 0.932. 
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Figure 8: The plot shows the nonmonotonic behavior of the particle density produced 
with the variation of the coupling constant. The value of b = M x /Hi is set to 1 here. 
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